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Abstract 

Concentrated solutions of short blunt-ended DNA duplexes, down to 6 base pairs, are known to 
order into the nematic liquid crystal phase. This self-assembly is due to the stacking interactions 
between the duplex terminals that promotes their aggregation into poly-disperse chains with a 
significant persistence length. Experiments show that liquid crystals phases form above a critical 
volume fraction depending on the duplex length. We introduce and investigate via numerical sim- 
ulations, a coarse-grained model of DNA double-helical duplexes. Each duplex is represented as an 
hard quasi-cylinder whose bases are decorated with two identical reactive sites. The stacking inter- 
action between terminal sites is modeled via a short-range square-well potential. We compare the 
numerical results with predictions based on a free energy functional and find satisfactory quanti- 
tative matching of the isotropic-nematic phase boundary and of the system structure. Comparison 
of numerical and theoretical results with experimental findings confirm that the DNA duplexes 
self-assembly can be properly modeled via equilibrium polymerization of cylindrical particles and 
enables us to estimate the stacking energy. 



PACS numbers: 64.70.mf,61.30.Cz,64.75.Yz,87.15.A-,82.35.Pq,87.14.gk 



I. INTRODUCTION 



Self-assembly is the spontaneous organization of matter into reversibly-bound aggregates. 
In contrast to chemical synthesis where molecular complexity is achieved through cova- 
lent bonds, in self-assembled structures the molecules or supramolecular aggregates sponta- 
neously form following the minimization of their free energy. Self-assembly is ubiquitous in 
nature and can involve the structuring of elementary building blocks of various sizes, ranging 
from simple molecules (e.g. surfactants) to the mesoscopic units (e.g. colloidal particles), 
thus being of topical interest in several fields, including soft matter and biophysics [TH3]. 
Understanding and thus controlling the processes of self-assembly is important in material 
science and technology for devicing new materials whose physical properties are controlled 
by tuning the interactions of the assembled components [H-0 171412] . 

A particular but very interesting case of self-assembly occurs when the anisotropy of at- 
tractive interactions between the monomers favors the formation of linear or filamentous ag- 
gregates, i.e. linear semi-flexible, flexible or rigid chains. A longstanding example is provided 
by the formation of worm-like micelles of amphiphilic molecules in water or microemulsions 
of water and oil which are stabilized by amphiphilic molecules . If supramolecular aggre- 
gates possess a sufficient rigidity the system may exhibit liquid crystal (LC) ordering even 
if the self-assembling components do not have the required shape anisotropy to guarantee 
the formation of nematic phases. An intense experimental activity has been dedicated to 
the study of nematic transitions in micellar systems [T3HTB] . Another prominent case is that 
of formation of fibers and fibrils of peptides and proteins [T6HT9] . Over last 50 years LC 
phases have been also observed in solutions of long duplex B-form DNA composed of 10 2 
to 10 6 base pairs [20H23] and in the analogous case of filamentous viruses [23H28]. More 
recently, a series of experiments [29H3T] ,have provided evidence that also a solution of short 
DNA duplexes (DNAD), 6 to 20 base pairs in length can form liquid crystal ordering above 
a critical concentration, giving rise to nematic and columnar LC phases [29|. 

This behavior was found when the terminals of the duplexes interact attractively. This 
condition is verified either when duplexes terminate bluntly, as in the case of fully comple- 
mentary strands shown in Fig. [Tk, or when the strands arrange in shifted double-helices 
whose overhangs are mutually interacting. This behavior is not restricted to B-form DNA 
oligomers, as it has also been observed in solutions of blunt-ended A-form RNA oligomeric 
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duplexes [32]. As terminals are modified to disrupt attraction, the LC long range order- 
ing is lost. Overall, the whole body of experimental evidence supports the notion that LC 
formation is due to the formation of reversible linear aggregates of duplexes, in turn pro- 
moting the onset of long-ranged LC orientational ordering. According to this picture, the 
LC ordering of oligomeric DNA is analogous to the LC ordering of chromonic liquid crystals 
[33] . Both in chromonics and in blunt-ended DNA duplexes, the aggregation takes place 
because of stacking interaction, generally understood as hydrophobic forces acting between 
the flat hydrocarbon surfaces provided by the core of chromonic molecules and by the paired 
nucleobases at the duplex terminals [3U ES] • 

The LC ordering of nucleic acids is relevant for various reasons. Firstly, it provides a 
new model of reversible aggregation leading to macroscopic ordering in which the strength 
of the inter-monomer attraction can be modified by changing the duplex terminals (bunt- 
end stacking or pairing of overhangs). Second, it provides a new access to the DNA-DNA 
interactions, and in particular to stacking interactions, whose nature is still investigated 
and debated [MJ ES] • In this vein, self-assembly acts as an amplifier of the inter-monomeric 
interactions, enabling studying the effects of minor molecular modification (e.g. oligomer ter- 
minations) on the base stacking. Finally, stacking and self-assembly are often invoked as the 
prebiotic route to explain the gap between the random synthesis of elementary carbon-based 
molecules and the first complex molecules, possibly RNA oligomers, capable of catalyze their 
own synthesis [US]- To proceed in any of these directions, it is necessary to rely on models 
enabling to quantitatively connect the collective behavior of nucleic acids oligomers to the 
molecular properties, and in particular to the duplex size and to the strength and range of 
the interduplex attractions. 

While the nematization transition in rigid and semiflexible polymers has been investi- 
gated in details in the past and rather accurate thermodynamic descriptions have been 
proposed [3"THlo] . much less is known for the case in which the nematic transition takes 
place in equilibrium polymer systems, i.e. when the average length of the chains depends 
on the state point explored. Recent theoretical and numerical works [SJ H7] has renewed 
the interest in this topic [H]. Ref. [UJ investigate the self-assembly and nematization of 
spheres, while Ref. [16] focuses on polymerization of interacting cylinders. In this article 
we propose a coarse-grained model similar to the one introduced in [IB] and devised to 
capture the essential physical features of equilibrium polymerization of DNA duplexes and 
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study it numerically via Monte Carlo simulations in the constant temperature and pressure 
ensembles, applying special biasing technique [4Uj I5T)] to speed up the equilibration process. 
We then develop a free-energy functional, building on Wertheim and Onsager [M] 

theories which provides, a satisfactory description of the system in the isotropic and ne- 
matic phases. A comparison of the calculated phase boundaries for different aspect ratio 
and different interaction strength with the experimental results allow us to confirm that 
the DNAD aggregation and LC ordering processes can be properly modeled via equilibrium 
polymerization of cylindrical particles and to provide an estimate of the stacking energy. 
In Section [IT] we introduce the coarse-grained model of DNADs and we provide some 



details of the computer simulations we performed. Section |III| gives a summary of the 
analytic theory which we developed in order to describe the system in the isotropic and 
nematic phases. A comparison of our analytical approach with numerical results is presented 



in Section [V] while in Section VI comparing our theoretical results with experimental data 



we provide an estimate of the stacking energy. In Section VII we draw the conclusions of 
our work. 




X=1 X=2 



FIG. 1: Coarse-grained model of DNA duplexes, (a) DNA duplex and a 3D graphical representation 
of its corresponding coarse-grained model comprising a SQ, symmetric around the x axis, decorated 
with two sticky spots located on its bases. The figure also show SQs of different aspect ratios 
(Xq = 1, 2, 3) and the projection of their surfaces onto the xy-plane. Note that the base roundness 
increases on increasing Xq. (b) A random chain of 10 monomers and a representation (blue 
clouds) of the points where the center of mass of a different monomer can be located in a bonding 
configuration. This set of points defines the bonding volume. 
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II. MODEL AND NUMERICAL DETAILS 



In this section we introduce a coarse-grained model devised to capture the essential phys- 
ical features of end-to-end stacking (equilibrium polymerization) of DNA duplexes and well 
suited to being investigated both theoretically and numerically In the model, particles 
(DNADs) are assimilated to superquadrics (SQ) with a quasi-cylindrical shape decorated 
with two reactive sites on their bases determining their interactions. SQs are a straightfor- 
ward generalization of hard ellipsoids (HE), their surface is in fact defined as follows: 
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where the parameters p, m, n are real numbers and a, b, c are the SQ semi-axes. In our case 
we set m = n = 2, p = 16 and b = c, so that the SQ resembles a cylinder with rounded 
edges (see Fig. [T]). Such SQs can be fully characterized by the elongation X = a/b and by 
the parameter p, that determines the sharpness of the edges (see Fig. [I]). SQs of elongation 
X < 1 are called "oblate", while SQs of elongation X > 1 are called "prolate", as for the 
HEs. As unit of length in our simulations we use the length of short semi-axes b. In the 
present study we investigated only prolate SQs with elongations X = 1, 2 and 3. We chose 
such elongations because DNADs used in experiments (29] have a diameter D = 2nm and 
are composed of 6 to 20 base pairs (BP) each 0.3 nm long, hence their elongation Xq ranges 
from 1 to 3. 

Each particle is decorated with two attractive sites, located along the symmetry axis 
(x-axis in Fig. [TJ at a distance d/b = X — 0.46 from the DNAD center of mass, in order to 
model hydrophobic (stacking) forces between DNADs. Sites belonging to distinct particles 
interact via the following square-well (SW) potential: 



Pusw = < ^ , (2) 




where r is the distance between the interacting sites, 5/b = 1.22 is the range of interaction 
(i.e. the diameter of the attractive sites), /3 = l/kgT and /eg is the Boltzmann constant. 
Therefore in the present model the anisotropic hard-core interaction is complemented with 
an anisotropic attractive potential in a fashion similar to what has been done in the past 
for other systems, like water [55J, silica [56J and stepwise polymerization of bifunctional 
diglycidyl ether of bisphenol A with pentafunctional diethylenetriamine [57J 158] . 



The location and diameter of the attractive sites have been chosen to best mimic the 
stacking interactions between blunt-ended DNAD and in particular they ensure that: 

1. the maximum interaction range between two DNADs bases is of the order of typical 
range for hydrophobic interactions (i.e. 2 A see [59J), i.e. of the order of water molecule 
dimensions 

2. the extent of the attractive surface of the DNADs bases is compatible with the surface 
of aromatic groups present in DNADs and which are responsible for hydrophobic 
interactions 

We note that in the present model each DNAD is symmetric around the x— axis (see Fig. 
[T|, hence we are neglecting rotations around it. 

We performed Monte Carlo (MC) simulations in the canonical and isobaric ensembles. 
We implemented the aggregation biased MC technique (AVBMC) developed by Chen and 
Siepmann |l9j [50] in order to speedup the formation of linear aggregates. To detect the 
overlap of two DNADs we calculated the distance using the algorithm described in Ref. 
[60] . In all simulations we adopted periodic boundary conditions in a cubic simulation box. 
We studied a system of iV = 1000 particles in a wide range of volume fractions and pressure 
P respectively. Initially we prepared configurations at high temperature with all DNADs 
not bonded, then we quenched the system to the final temperature (i.e. to the final value of 
(3AEs) letting the system equilibrate. We checked equilibration by inspecting the behavior 



of the potential energy and the nematic order parameter (see Sec. VB) in the system 



III. THEORY 



Following the work of van der Schoot and Cates [HJ HE] and its extension to higher volume 
fractions with the use of Parsons-Lee approximation [611 [62] as suggested by Kuriabova et 



6 



al. we assume the following expression for the free energy of our system: 



v 



oo 



J>(0 {In KKO] - 1} + 



1=1 



+ 



imj2v(i)v(i')v exd (i,n 



1=1 
l'=l 



oo 



oo 



(PAE S + a h ) - l)v{l) + "(l)°o(l) 



(3) 



i=i i=i 



where v(l) is the number density of chain of length I, normalized such that Ylk=i I KO = Pi 
Vd is the volume of a monomer, (3AEs is the (positive) stacking energy, v exc i(l,l') is the 
excluded volume of two chains of length I and I' and <Jb is the entropic free energy penalty 
for bonding (i.e. is the contribution to free energy due to the entropy which is lost by 
forming a single bond) and r\ (</>) is the Parsons-Lee factor [6T] 



and <r [l5] accounts for the orientational entropy that a chain of length I loses in the nematic 
phase (including possible contribution due to its flexibility). Differently from Ref. [46, 48J, 
but as in Ref. [17], we explicitly account for the polydispersity inherent in the equilibrium 
polymerization using a discrete chain length distribution. We explicitly separate the bonding 
free energy in an energetic (/3AEs) and entropic (er&) contribution. Differently from Ref. [471 
14*8] . but as in Ref. [46] we include the Parson-Lee factor. Indeed, the Parson decoupling 
approximation satisfactory models the phase diagram of uniaxial hard ellipsoids [63], hard 
cylinders j5"4"] . linear fused hard spheres chains [(53], mixtures of hard platelets [SS], hard 
sphero-cylinders [5TH5T'] . rod-plate mixtures [70], mixtures of rod-like particles [TTJ [72] and 
mixtures of hard rods and hard spheres [73]. On the other hand, Ref. [74] finds that the 
Parsons theory is not satisfactory in the case of rigid linear chains of spheres. 

A justification of the use of Parsons-Lee factor in Eq. ^ for the present case of aggre- 
gating cylinders is provided in Appendix A.Here we only note that the present system, in 
the limit of high T where polymerization is not effective, reduces to a fluid of hard quasi- 
cylinders, where the use of Parsons-Lee factor is justifiable[641 [68] 169] . Moreover, in the 
dilute limit (r/(</>) — > 1) the excluded volume term in Eq. [3] reduces to the excluded volume 
of a polydisperse set of cluster with length distribution which is conform to Onsager 
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original theory [54]. In other words the form chosen in Eq. [3] for the excluded volume contri- 
bution to the free energy reduces to the correct expressions in the limit of high temperatures 
and in that of low volume fractions. 

Following Van der Schoot and Cates [HI SB] a generic form for v exc i(l, I') can be assumed 
as a second order polynomial in I and Z', with 



*w[U';/(u)] = 2 J /(u)/(u')£> 3 [*!(7,X )+ 



+ / -±-^ 2 ( 7 ,X )X + * 3 (7,*o)*o H']dfldfl' (5) 

where /(u) is the probability for a given monomer of having orientation u within the solid 
fi and O + dfl and ty a describe the angular dependence of the excluded volume. The 
orientational probability /(u) are normalized as 



+ -E(sin7))^l + 2X 2 sin7 1 1' 

7T 2 



f(n)dQ = 1 (6) 

In particular for two rigid chains of length I and V which are composed of hard cylinders 
(HC) of diameter D and length X D v exc [(l, V) has been calculated by Onsager in 1949 

<W(M') = J f(u)f(u')D 3 [| s in 7 +|Xo(l + |cos7| + 

dSldQ! (7) 

where cos 7 = u ■ u' and E{sm r )) is the complete elliptical integral 

1 C 2n 

E(sinj) = - (1 - sin 2 7 sin 2 ^) 1/2 # (8) 
4 Jo 

On passing we observe that the integrals in Eq. ([7| can be calculated exactly in the isotropic 
phase while in the nematic phase the calculation can be done analytically only with suitable 
choices of the angular distribution /(u). Comparing Eqs. ([7]) and ^ for HC one has: 

7T 

^1(7,^0) = 2 sin ^ 

#2(7, *o) = + |cos 7 | + -E(sin 7 )) 

Z 7T 

^3(7, X ) = 2 sin 7 (9) 

In view of Eqs. ^ we notice that for HCs the functions ^1(7), ^2(7) and ^3(7) accounts for 
the orientational dependence of the excluded volume of two monomers having orientations 



u and u' with u • u' = cos 7. It is also worth observing that the first term of the integrand in 
Eq. ([7]) is independent of I and hence accounts for the excluded volume interaction between 
two HCs end caps, the second term is linear in I and V and accounts for the excluded volume 
between the two end caps of a chain and all midsections of the other one and the third 
one proportional to IV models the interaction between all IV pairs of midsections of the two 
chains [TH HH]. In summary, Eq. (|5| is exact for two rigid chains of HCs but according to 
[El [39] to lowest order of approximation it is justifiable to use such equation also for two 
semi-flexible chains. We then assume that v exc i remains additive with respect to end-end, 
end-midsection and midsection-midsection excluded volume contributions even if the chain 
is semi-flexible. Finally our further ansatz is that Eq. ^ is also a good functional form for 
the excluded volume of two superquadrics having quasi-cylindrical shape: we will check the 
validity of this hypothesis using our simulations data. 

An exact expression for a is not available. The two following limits have been calculated 
by Khokhlov and Semenov SOI US] : 

< * w =i/(f) , '- ldn+ 
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Finally, we notice that in the limit of rigid rods with fi(u) = f(u)v(l), (the same limit 
selected in Ref. [46J), the free energy in Eq. (|3| reduces to: 
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which is analogous to the free energy expression used by Glaser et al. 



A. Isotropic phase 



In the isotropic phase all orientations are equiprobable, hence: 
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Plugging Eq. ( 12 ) into Eq. ^ and calculating the integrals one obtains: 
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For hard cylinders the excluded volume can be calculated explicitly and it turns out to 



be: 
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Building on Eq. (14), the generic expression for the excluded volume v exc i(l,l') reported 



in Eq. ^ in the isotropic phase takes the form: 

Vexcl(lJ') = 2 



Aj(X ) + k I (X )v/-^+ 



BjiX^lV 



(15) 



We assume that the chain length distribution v(l) is exponential with a average chain 
length M 

u(l)= P M~ {l+l \M -I) 1 ' 1 (16) 

where 
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With this choice for v(l) the free energy in Eq. (13) becomes: 
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Note that in general ki, Bj and Ai depend on Xq. 

The minimum of the free energy with respect to M (i.e. the equation d(/3F/V)/dM = 
provides the searched equilibrium value for M. Dropping terms in 0(1/M 2 ) one obtains 

1 
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(19) 



where u = 4e (Jb . This formula for M differs from the one reported by Kindt [17] by the 
presence of the Parsons-Lee factor, which will play a role at high volume fractions. 



The expression for M in Eq. ( 19 ) coincides with the parameter-free expression for average 



chain length M w obtained within Wertheim's theory (e.g. see Refs. [5TH531 1751 176]). when 
is small and e kl ^ v ^ ~ 1. Indeed, in Wertheim theory 



U„ =I + Iwi + 8^A 
2 2V v d 
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where A = Vfe(e /3A£ ' S — 1) and Vf, is the bonding volume [7S]. In the limit e^ AEs 3> 1, always 
valid in the T-region where chaining takes place, 



M w = \ + \Jl + %%e^ E s 
2 2V Vd 

The equivalence between the two expressions provide an exact definition of u> as 



(21) 
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Although Eq. (19) has been derived ignoring 0(1/ M ) terms in the free energy, the 



average chain length M can be always calculated, and this is what we do in this work, 
finding numerically the zero of d((3F/V)/dM = 0. 
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B. Nematic Phase 



In the nematic phase the function /(u) depends explicitly on the angle between a given 
particle direction and the nematic axis, i.e. on the axis u. The function /(u) is also 
called the "trial function" and it generally depends on a set of parameters that have to be 
obtain through the minimization of the free energy. Also in the nematic phase we assume 
an exponential distribution for 2/(7). In view of the analytical expression for the excluded 
volume v ex d for cylinders, we assume the following form for the v exc i of two DNADs averaged 
over the solid angle using a one parameter (a) dependent trial function 



v exd {l,l' ,oc) 
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(23) 



If we insert Eqs. (57), (23) and (16) into Eq. ^ we obtain after some algebra: 
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where a Q = Y^i a o{l)v(l) 



C. Phase Coexistence 



P 2 + 



(24) 



Using the free energy functionals in Eqs. (18) and (24) the phase boundaries, i.e. 0^ = 
VdPN an d <pi = VdPii of isotropic- nematic transition can straightforwardly calculated by 
minimizing the free energy with respect to average chain lengths in the isotropic and nematic 
phases, i.e. Mi and Mat, and a. We also require that the isotropic and nematic phases have 
the same pressure, i.e. Pi = Pn and the same chemical potential \ii = //jv- These conditions 
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require numerically solving the following set of equations: 

d -F iS0 ( Pl , Mj) = 



dMi 
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Kem(PN,M N ,a) = 



dM N 
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—F iso { Pl ,M h a) = 
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P/(p/,M 7 ) = P N (p N ,M N ,a) 

(j,i(pi, Mi) = fi N (p N ,M N ,a) (25) 

IV. CALCULATION OF FREE ENERGY PARAMETERS 

The theory illustrated in the previous section requires the calculation of several parame- 
ters, Vb, kj, A T , Bf, fc/v, A N , B N , l p . Since for super-quadrics an explicit calculation of these 
parameters is very unlikely in the following we describe simple methods to calculate them 
numerically. For example the calculation of the excluded volume between clusters and the 
calculation of the bonding volume require the evaluation of complicated integrals, which can 
be estimated with a Monte Carlo method. The general idea indeed behind Monte Carlo is 
that such complicated integrals can be calculated by generating a suitable distribution of 
points in the domain of integration. 

A. Excluded volume in the isotropic phase 



In the isotropic phase, v exc i(l, V) can be written as reported in Eq. (15). If I = I', 

v exd (l,l)=2Aj + 2k lVd I + 2BjX 2 I 2 (26) 

Hence, from a numerical evaluation of v exc i(l, I) for several I values (whose detailed procedure 
is described in Appendix [x| it is possible to estimate Aj, kj and Bj. Fig. [2] shows v exc i(l, I) /I 
vs I. A straight line describe properly the data for all Xq values, suggesting that Aj 0. 
From a linear fit one obtain 2BjX^ (slope) and 2kjVd (intercept). 

B. Calculation of the Bonding Volume 

The bonding volume Vb can be calculated numerically performing a Monte Carlo calcu- 
lation of 
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FIG. 2: (a) Excluded volume of two chains of length I calculated numerically as a function of I for 



Xq = 1, 2, 3. Dashed lines are fits to Eq. (26). (b) Bonding volume as a function of elongation Xq. 



V b = J 9 (-AE S - u sw - V HC ) dv dtoi £l 2 (27) 

where Vhc — VHc{ r , S^i, ^2) is the hard core part of the interaction potential and 8(x) 
is the Heaviside step function, i.e. 9(x) = 1 if x > or otherwise. The detail of the 
numerical integration are reported in Appendix |X} The resulting values of V& for different 
X are shown in Fig. [2^b). V& grows with X , an effect introduced by the different rounding 
of the SQ surface close to the bases. Indeed, how is it shown in Fig. [TJ on increasing Xq the 
base surface is more rounded and such different rounding offers a different angular width 
over which bonds can form, an effect which will also reflect in the Xq dependence of the the 



persistence length of the self-assembled chains, as it will be discussed in details in Sec. IV E 
The dependence of the bonding angle on Xq is not present for HCs and in that case the 
bonding volume would be constant. 
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C. Nematic phase 



In what follows we assume that the angular distribution of the particle orientation can 
be well described by the one parameter Onsager trial function[54j, i.e.: 



a 

f(u) = f (u) = — cosh(acosg) (28) 
47r smh a 

where 9 is the angle between the particle and the nematic axis and the system is supposed to 
have azimuthal symmetry around such axis. The excluded volume v exc i(l,l,a) between two 
clusters of equal length I can be calculated using the procedure illustrated previously for the 
isotropic case with the only difference that now monomers are inserted with an orientation 



extracted from the Onsager angular distribution defined in Eq. (28). 



To estimate numerically An (a), &tv(«) and Bn{ol) we specialize Eq. (23) to the case of 
I — I' — 2, I — I' — 3, I — I' — 4, evaluating numerically for several values of a v exc i(2, 2, a), 



fW«(3,3,a) and v exc i(4 } 4, a). Inverting Eq. (23) allows us to express A N (a), k^ia) and 
B N (a) as a function of v exc i{2, 2, a), v exc i(3, 3, a) and v exc i(4:, 4, a) as explained in detail in 
Appendix C. 



D. Estimate of the orientational entropy in the nematic phase 

We propose to model the orientational entropy in the nematic phase using the following 
expression proposed by Odijk jlS] (other possibilities can be found in Refs. [77j and [T8"] ) 
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This expression, in the limit of "rigid chains" (RC) and "flexible chains" (FC) reduces 
approximately to the exact limits [39J 
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Latter formulas can be also obtained plugging the Onsager trial function /o( u ) in Eqs- 
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Unfortunately, Eq. (29) is hardly tractable in the minimization procedure requested to 



evaluated the equilibrium free energy and hence the two approximations in Eq. (30) are 
often preferred. While in the case of fixed length polymers, the knowledge of the persistence 
length selects one of the two expressions, in the case of equilibrium polymers, different chain 
length will contribute differently to the orientational entropy. In particular, when the chain 
length distribution is rather wide, it is difficult to assess if the RC (chosen in Ref. jl6]) or 
the FC (chosen in Ref. [J7]) limits should be used. To overcome the numerical problem, 
still retaining both the RC and the FC behaviors, we use the following expression for a : 

l=l -l s 

to = J2 "(0 - x i + 

1=1 ^ 

I — Iq 

in which the contribution of chains of size Iq is treated with the RC while the contribution of 
longer chains enters with the FC expression. We pick Iq by requesting maximum likeihood 



between Eq. (31) and Eq. (29) in the relevant M-a domain. We found that the value Iq ~ 9 



is appropriate for most studied cases. 



E. Estimate of persistence length 



In order to estimate the persistence length, entering in Eq. (31 ), we randomly build chains 
according to the procedure described in Appendix BWe estimate the "chain persistence 
length" l p by evaluating following spatial correlation function: 

C (\i-j\)^J2(±(i)-m) (32) 

hi 

where i,j label two monomers along the chain (i = is the first monomer at chain end) 
and x(z) is a unit versor directed along x-axis of the monomer (i.e. their axis of symmetry, 
see Fig. [T]), that coincides with the direction along which the two attractive sites lie. (...) 
denotes an average over the whole set of independent chains which has been generated. 

In Figure|3]we plot Co(\i — j\) for all elongations studied. All correlations decay following 
an exponential law, whose characteristic scale is identified as the persistence length (in unit 
of monomer). In the explored Xq range, 10 < l p < 25. The more elongated monomers have 
a smaller persistence length. The X dependence of l p arises from the different roundness of 

16 



the bases (implicit in the use of SQ), as discussed in the context of the bonding volume and 
in Fig. [Tj 
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FIG. 3: Spatial correlation function Co(\i — j\) (see text for its definition) calculated generating 
random chains of 50 monomers for aspect ratios Xq = 1, 1.5, 2, 2.5, 3. Dashed lines are fits to the 
functional form Co(\i — j\) = exp[— \i — j \ /l p ]. From these fits the chain persistence length l p can 
be estimated (see legend). 



V. RESULTS AND DISCUSSION 



In this section we compare results from simulations with theoretical calculations based 



on the theory discussed in Sec. Ill 



A. Isotropic phase 



Fig. [4] (a)-(c) show the packing fraction dependence of M for Xq = 1,2,3 for all tem- 
peratures investigated. The dashed curves are calculated by minimizing with respect to M 



the isotropic free energy in Eq. (18) using the values of V&, ki and Bj obtained in Sec. 
IV A without any fitting parameter. Up to volume fractions around ~ 0.20 the agreement 
between theoretical and numerical results is quite good for all cases considered. Above such 
volume fraction the theoretical predictions start deviating appreciably, a discrepancy that 



17 



we attribute at moderate and high to the inaccuracy of the Parsons decoupling approxi- 
mation. 

In Fig. p] (d) we report the cluster size distribution v(l) as obtained from both simulation 



and theory, the latter calculated according to Eq. ( 18 ) with M obtained by minimization of 
the isotropic free energy. As expected, the cluster size in the isotropic phase is exponential. 
These results suggest that a reasonable first principle description of the isotropic phase is 



provided by the free energy of Eq. (18), when the parameter of the model are properly 
evaluated. 



B. Nematic phase 

On increasing the system transform into a LC phase. We estimate the degree of nematic 
ordering by evaluating the largest eigenvalue S of the order tensor Q, whose components 
are: 

1 x 3 1 

Qa/3 = ~jy -((Ui) a ( Ui )^) - -5 aS (33) 

i 

where a(3 G {x,y,z}, and the unit vector (uj(t)) a is the component a of the orientation 
(i.e. the symmetry axis) of particle i at time t. A non-zero value of S signals the presence 
of orientational order in the system and it can be found not only in the nematic phase but 
also in partially ordered phases as columnar and smectic phases. Since in this article we 
focus only on the nematic phase, to verify that the simulated state points are not partially 
ordered we calculate, following Ref. jl6], the three dimensional pair distribution function 
g(r) defined as: 

f( r ) = ^(EE*( r -( r *- r i))) ( 34 ) 

\ i=l j^i I 

where 5(x) is the Dirac delta function. We calculate the g(r) in a reference system with 
the z-axis parallel to the nematic director. Figure [5] shows g(x,y,0) and g(0, y, z), which 
correspond respectively to the correlations in a plane perpendicular to the nematic director 
and in a plane containing it for a given nematic state point ( Xo = 2, <fi = 0.38, (3AEs = 
8.33)). The g(x,y,0) is found to be isotropic, ruling out the possibility of a columnar 
or crystal phase (no hexagonal symmetry is indeed present). The g(0,y,z) reflects the 
orientational ordering along the nematic direction and rules out the possibility of a smectic 
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FIG. 4: (a)-(c) Average chain length M against <f> for Xq = 1,2,3 at all values of f3AE$ studied. 



Dashed lines are theoretical predictions calculated by minimizing the free energy in Eq. (18) with 



values of Vj,, kj and Bj derived by the procedure described in Sec. IV (d) Cluster size distributions 



(colored symbols) for several state points together with theoretical predictions (dashed lines). 



phase (no aligned sequence of peaks are present [IE]). Fig. |5|c) shows also a snapshot of the 
simulated system at the same state point. 

In what follows, we have systematically calculated and inspected g(r) to verify that all 
state points having a value of S large enough to be considered nematic are indeed transla- 
tionally isotropic, i.e. with no translational order. 

Fig. [6] shows the nematic order parameter and the average chain length M calculated 
from simulations as well as with the theoretical methodology described previously for two 
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different elongations at PAEs = 8.33. The theoretical value for S is obtained according to: 

/"X pns 2 f) — 1 
2ti- — f (9; a) sine d0 (35) 

Fig. [6] (a) shows that the nematic order parameter is very well captured by the theory, 
while the average chain length shows a clear disagreement between theory and simulations, 
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FIG. 6: (a) Average chain length and nematic order parameter S for several nematic state points. 
Dashed lines with stars (theory IM) are theoretical predictions assuming that monomers are 
isotropic (see text for details), (b) Cluster size distribution for two state points (Xq = 2, <p = 
0.38, PAE S = 8.33) and (X = 3,0 = 0.34, pAE s = 8.33). Circles are numerical results and 
dashed lines are exponential fits. The inset shows the chain length dependent nematic order pa- 
rameter Si for the same state points. 



again suggesting that the error introduced by the Parsons decoupling approximation which 
was already observed in the isotropic case for large packing fractions is here enhanced by the 
further increase in cf). Another possible source of error could arise from the hypothesis that 
the cluster size distribution is exponential also in the nematic phase. To test this hypothesis 
we show in Fig. [6] (b) the cluster size distributions in two different state points. In all cases, 
the distributions are not a single exponential. This phenomenon has been already observed 
and discussed by Lu and Kindt jl?] and described as a two exponential decay of v(l) with 
the exponential decay of short chains extending up to I ~ 50. They took into account such 
bi-exponential nature of the distribution to better reproduce the isotropic-nematic phase 
boundaries |79j in their theoretical approach. In the present case, only very short chains 
(not to say only the monomers), fall out of the single exponential decay. To test if the 
different decay reflects a different orientational ordering of the small clusters compared to 
long chains, we follow Ref. [79J and evaluate the length-dependent nematic order parameter 
Si, that is the nematic order parameter calculated for each population of clusters of size 
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I. The results, reported in the inset of Fig. |6Jb), show that Si is around 0.7 — 0.8 for all 
clusters sizes except for I = 1, i.e. except for monomers. 

To assess how much the theoretical predictions are affected by the assumption of a single 
exponential decay (and of the associated identity of S for all chains), we evaluate the cor- 



rection of the free energy functional in Eq. (24) arising from the assumption that monomers 



are isotropic, while all other chains are nematic. To do so we exclude from the calculation 
of the orientational entropy the monomers and take into account in the calculation of the 
excluded volume contribution the fact that monomers are isotropic. The revised free energy 
can be thus written as 
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Minimizing such expression we calculate the resulting improved estimate for the average 
chain length and nematic order parameter and the results are also shown in Fig. |6j The 
new estimates slightly improve over the previous ones, suggesting once more that the leading 
source of error in the present approach, as well in all previous, has to be found in the difficulty 
of properly handling the term successive to the second in the virial expansion. 
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C. Phase Coexistence 



NPT-MC simulations provide a rough estimate of the location of phase boundaries, being 
affected by the hysteresis associated to the metastability of the coexisting phases. It is 
possible thus only to bracket the region of coexistence, by selecting the first isotropic state 
point on expansion runs which started from a nematic configuration and the first nematic 
state point on compression runs started from an isotropic configuration. We performed 
NPT-MC simulations for X = 2 and (3AE$ = 6.67 in a wide range of pressures P for a 
system of 1000 SQs. The resulting equation of state is shown in Fig. [7ta). As expected 
a clear hysteresis is observed, which allows us to detect some overestimated boundaries for 
the isotropic-nematic transition. The same figure also reports the theoretical estimates of 
the transition. The theoretical critical pressure is smaller than the numerical one, resulting 
into a extended region of coexistence then numerically observed. Comparing the values of 
the pressure predicted by the theory with the simulation values, we notice that the main 
error arises from the pressure of the nematic phase which is underestimated. Finally, Fig. 
7] (b)-(d) show the predicted phase diagram for several values of (3AE$ as a function of the 
elongation. On increasing (5AEs (i.e. decreasing T or increasing the stacking energy) there 
is a small decrease of <pi and a significant decrease of 0jv, resulting in an overall decrease 
of the I — N coexistence region. Such trend can be understood in term of increase of the 
average chain length resulting from the increase of (3AE$- As expected, both <pi and 4>n 
decrease on increasing X . 

Finally we recall that in our model the persistence length l p depends on the elongation 



as discussed in Section IV E , anyway the dependence on l p of the theoretical phase diagrams 



showed in Fig. [7j (b)-(d) is negligible at the present level of accuracy of our theoretical 
calculations. 



VI. COMPARISON WITH EXPERIMENTS 

Refs. [29] and [31] report the critical concentrations (c), in mg/ml, for the I—N transition 
of blunt-ended DNAD. These experimental data can be transformed into volume fractions 
once the relevant properties of DNAD are known (DNAD molecular weight = 660 Aft, 
diameter D « 2 nm, length L = iVft/3 nm, where Nj, is the number of bases in the sequence). 
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FIG. 7: (a) Equation of state {P vs <f>) calculated compressing an isotropic initial configura- 
tion (squares) or expanding an initial nematic configurations (circles). Vertical dashed lines show 
the theoretical predictions for the phase boundaries, (b)-(d) Coexistence regions predicted from 
theoretical calculations for the 3 stacking energies values PAE$ investigated. 



The number density p of DNADs is related to the mass concentration 

P = TT— ( 39 ) 
F Nm D v ' 

Since = LD 2 ir/4: is the volume of a DNAD, the volume fraction can be expressed as: 

cLD 2 n . . 

^ = pVd = W^ (40) 

Data in Refs. [29] and [31] put in evidence that blunt-end duplexes of equal length but 
different sequences may have different transition concentrations. As discussed in Ref. [31] 
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this phenomenon can be attributed to the slight differences in B-DNA helical conformation 
resulting from the difference in sequences. These differences induce some curvature in the 
DNAD aggregates, in turn enhancing the transition concentration. Indeed, sequences that 
are known to form straight double helices order into the N phase at lower concentrations. 
Therefore, for each oligomer length in the range 8-16 bases, we selected the lowest transition 
concentration among the ones experimentally determined, since these would be relative to 
duplexes closest to the symmetric monomers considered in the model. Such values have been 
reported in Figure [8] as a function of base number ]V& (top axis) and as a function of X 
(bottom axis). Apart for iV& = 12, of which a large number of sequences have been studied, 
the transition concentrations for the other iVj, values would probably be corrected to lower 
values if a larger number of sequences were experimentally explored. We would expect this 
to be particularly true for the shortest sequences, in which the effect of bent helices could 
be more relevant. 

In the experiments, DNADs are in a water solution with counter-ions resulting from 
the dissociation of the ionic groups of the phosphate-sugar chain. Given the high DNA 
concentration necessary for the formation of the N phase, corresponding to concentration of 
nucleobases in the 1M range, the ionic strength simply provided by the natural counter-ions 
is large enough to effectively screen electrostatic interactions between DNADs. This becomes 
less true for the longest studied sequences, for which the transition concentration is lower. 
We hence decided to perform a small number of test experiments on the iV& = 20 oligomers 
with a double purpose: (i) determine more accurately the transition concentration value for 
this compound and (ii) test the effect of varying the ionic strength predicted by the model 
here described. With respect to a fully screened DNAD where electrostatic repulsion can be 
neglected, a partly screened DNAD has a larger effective volume, thus filling a larger volume 
fraction of the solution, and a smaller axial ratio X , since electrostatic repulsion is equal in 
all directions. Therefore, adding salt would bring about two competing effects: the reduction 
in particle volume enhances the concentration needed to reach the I — N phase boundary, 
while the grown of Xq could favor the nematic ordering even at lower concentrations. 



In particular according to Eq. (40) the following relation between the critical concentra- 
tion cn and the critical volume fraction (j)^ holds: 

c N = MXo)^ (41) 
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On the basis of the phase diagrams of Figure^ (b)-(d) 4>n(Xq) depends weakly on X = L/D, 
i.e. 4>n(X ) w (p N) where (j>° N is a constant. Hence the theory introduced in the present paper 
predicts that a reduction of DNAD effective volume due to the addition of salt (i.e. a decrease 



of LD in Eq. (41 ) leads to an overall increase of the concentration required for N ordering. 

We have measured the transition concentration of the self-complementary 20mer 
CGC GAAAATTTTC GC G , a sequence whose I — N transition at room temperature was 
previously measured and determined to be cjn ~ 200 rag j ml [29]. With the same method, 
based on the measurement of the refractive index of the solution, we determined the I — N 
transition concentration at room temperature at three different ionic strengths. The val- 
ues we obtained are c^v ~ 215mg/ml (no added salt), c^v ~ S20mg/ml (0.8 M NaCl), 
c/7v ~ 380mg/ml (1.2 M NaCl). The data indicate that the onset of the nematic order- 
ing in solutions of 20mers is indeed sensitive to the ionic strength, and that the transition 
concentration grows upon increasing the amount of salt, as expected on the basis of our 
theoretical calculations for the present model. In Figure [8] we display the transition volume 
fraction derived by the transition concentration measured for 1.2 M NaCl. At this ionic 
strength, the total concentration of Na + (dissociated from the oligomers + added with the 
salt) is about the same as the one resulting from counterions dissociated oligomers in the 
more concentrated solutions of shorter (8-12 mers) oligomers. 

Figure [8] compares the experimentally determined transition volume fractions with the 
values calculated from the model for (3AEs = 6.67 and f3AEs = 5.56. Although experi- 
mental data are noisy, they fall in the range AE$ ~ 5 — 7 (in units of kgT). Despite all 
the simplifying assumptions and despite the experimental uncertainty, results in Figure [8] 
provide a reasonable description of the X dependence of 4>n- 

In comparing the model with the experimental results, it is necessary to take note of the 
fact that the stacking energy between nucleobases, and thus the interaction energy AE S 
between DNAD, is temperature dependent, i.e. its entropic component is relevant [35J. 
This is a general property of solvation energies and thus it is in line with the notion that 
stacking forces are mainly of hydrophobic nature. Therefore, the range of values for AE$ 
determined in Figure [8] should be compared to the values of AG for the stacking interactions 
at the temperature at which the experiments were performed. Overall, the estimate of AE$ 
here obtained appears as in reasonable agreement with the free energies involved in the 
thermodynamic stability of the DNA double helices and confirms the rough estimate that 
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was given before (see the supporting online material associated to Ref. [29J ) . 
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FIG. 8: Critical volume fractions $jy as a function of elongation Xq (or equivalently iVj,) from 
theoretical calculations (for f3AE$ = 6.67 and f3AE$ = 5.56) and experiments [29J (circles). 



VII. CONCLUSIONS 

In this study we have developed a free energy functional in order to calculate the phase 
diagram of bi-functional quasi-cylindrical monomers, aggregating into equilibrium chains, 
with respect to isotropic-nematic transition. The model has been inspired by experiments 
on the aggregation of short DNA which exhibit at sufficiently high concentration nematic 
phases and the comparison between the theoretical predictions and the experimental results 
allows us to provide an estimate of the stacking energy, consistent with previous propositions. 

Our approach is quite general, parameter free and not restricted to particular shapes. We 
provide techniques to evaluate the bonding volume and the excluded volume, which enters 
into our formalism via the Parson-Lee decoupling approximation. We build on previous work, 
retaining the discrete cluster size description of Ref. [47J and the Parson-Lee factor for the 
excluded volume contribution proposed in Ref. [IB]. With respect to previous approaches 
we (i) explicitly account for the entropic and energetic contributions associated to bond 
formation; (ii) we do not retain any adjustable fit parameter. 

The resulting description of the isotropic phase is rather satisfactory and quantitative up 
to cf) ~ 0.2. The description of the nematic phase partially suffers from some of the approx- 
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imations made in deriving the free energy functional. More specifically, several signatures 
point toward the failure of the Parsons decoupling approximation in the (f) range typical of 
the nematic phase. While there is a sufficient understanding of the quality of such approxi- 
mation for monodisperse object (63J, EU |66j EHl EHJ [73], work need to be done to assess the 
origin of the failure of this approximation in the equilibrium polymer case and to propose 
improvements. 

We finally remind that the model here introduced does not consider the azimuthal rota- 
tions of each monomer around its axis. This neglect is adequate when the aggregation does 
not entail constraints in the azimuthal freedom of the monomers. This is the case of base 
stacking, in which the angular dependence of the stacking energy is arguably rather small. 
However, this is not the case of DNAD interacting through the pairing of overhangs and 
of the LC ordering of RNA duplexes. Because of its A-DNA-type structure, the terminal 
paired bases of RNA duplexes are significantly tilted with respect to the duplex axis, thus 
establishing even in the case of blunt-ended duplexes a link between the azimuthal angle 
of the aggregating duplexes and the straightness of the aggregate. However, with minor 
modifications the model here introduced could become suitable to include these additional 
situations. The limiting factor in developing such extension is the lack of knowledge to 
quantify the azimuthal constraints implied by these interactions. This situation, as well as 
the effects of off-axis components of the end-to-end interduplex interactions, will be explored 
in a future work. 
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IX. APPENDIX A 

Here we provide a justification for the use of Parsons decoupling approximation in the 
case of linear chains poly-disperse in length (with distribution based on the extension 

of Onsager's second- virial theory to mixtures of non-spherical hard bodies proposed in Ref. 
[80] . The contribution F exc i to the free energy due to excluded volume interactions between 
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chains can be written if we neglect intrachain interactions [80j EE] : 



V 6 J -J -J V ^ 

^( r , ni, n 2 )/(no/(n 2 ) r • v r F HC (r, n a , n 2 ) (42) 

where r is the distance between the centers of mass of the two chains 1 and 2, ill = 
{uj, . . . u[} and f2 2 = { u 1, • • • u y} are the orientations of the two chains, where uf is the 
orientation of monomer % belonging to chain a = 1, 2, gw(r, fii, f2 2 ) is the molecular radial 
distribution function of the mixture, which represents the correlations between two chains 
of length I and I', whose relative distance is r and which have orientations fii and fi 2 
respectively, Vnc{ r , S^i, ^2) is the hard-core part of the interaction potential and f(tt a ) 



is the angular distribution function of chain a. We note that in Eq. (42) the integration 
in p' is performed keeping fixed all the parameters related to f(Q a ). Neglecting intra- 
chain interactions is equivalent to ignore self-overlaps of chains, an assumption which is 
appropriate if chain length is not much greater than its persistence length and the chains 
can be considered non-extensible. 

Parsons decoupling approximations in this case accounts to putting: 

g u ,(r, n u n 2 ) = g§ S [r/a w (f, fl u n 2 )] (43) 

where g§ s is the radial distribution function of a mixture of hard spheres and <Ju>(r, fix, fi 2 ) 
is an angle-dependent range parameter which depends on chain lengths I and V . If the pair 
interaction is of the special form 

V^c(r,n!,n 2 ) = V H c[r/<T lv (r,tox,to2)] (44) 



noting that r • V r = rj^, Eq. (42) becomes: 
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With the substitution y = rjuiy from Eq. (45) one obtains: 
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The derivative of Vhc is a delta function hence we need only to evaluate the value of 
9u' S (y) a ^ contact (i.e. y 



and Eq. (46) becomes: 
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This expression tends to Parson's expression when the system is monodisperse {y{£) = p5n 
In the specific case of spherical particles au>(r, fix, ^2) = o"(r, fix, fl 2 ) = a and 
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(47) 



i.e. the excluded volume of two spheres of diameter a. Hence we are allowed to make the 
identification: 



V e xcl{l,l') 



dr dflx dfl 2 -f(flx)f(fl 2 )4>(r, ^1, fi 2 ) 



and write: 
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(48) 



(49) 



We note that the identification made in Eq. (48) can be also further justified using the same 
reasonings given in Sec. Ill As discussed in Ref. (80] a possible expression for g^, s is the 



one derived by Boublfk [82], which generalizes the Carnahan-Starling relation [83] for pure 
hard spheres to the case of mixtures, i.e. 
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(50) 



1 — C3 ' (1 - C3) 2 On + o-fi' 
where an is the diameter of an hard sphere corresponding to a chain of length I and ( n = 
(7r/6) J2i v {l)®?i- To map the system of polydisperse chains onto the equivalent mixture of 
hard spheres we need an expression for an. According to Ref. [80], the simplest choice is to 
consider spheres having the same volume of the corresponding linear chain of length I, i.e. 



v d 
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-ar 



(51) 



where we recall that Vd is the volume of a monomer. Although in principle we could use Eq. 



(49) together with Eqs. (50) and (51) to calculate the free energy contribution due to the 



excluded volume between particles, if we make the further assumption that 



HSi 



(52) 
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i.e. if we approximate the radial distribution function of the hard spheres mixture at contact 
with that of a monodisperse system of hard spheres having the same total volume fraction 
(i.e. setting in Eq. (50) aw = o with Mva = (n/6)a 3 ), we finally obtain 

§^ = m^ vmnvndM (53) 

IV 

where we used the Carnahan-Starling expression for g p') and we performed the in- 



tegration in p' . Eq. (53) is exactly the expression for the contribution to the free energy 



due to steric repulsion which we used in Section III In summary according to the above 



derivation we argue that Eq. (53) can be not accurate at high volume fractions due to the 



approximations made in Eqs. (43) (i.e. the Parsons decoupling approximation) and (52). 



Within the present treatment Eq. (53) is also not appropriate for chains with I ^> l p because, 
as already noted, chain self-overlaps can be significant and the hard body pair potential Vhc 



does not have the special form assumed in Eq. (44). 



We finally note that the approximation made in Eq. (52) can be avoided if one resorts to 



Eq. (49) instead of Eq. (53), although the required free energy calculations would become 



much more complicated. Anyway we verified for the isotropic phase that employing Eq. (49) 



instead of Eq. (52) does not provide any appreciable improvement in the present case. 



X. APPENDIX B 



The procedure to calculate the excluded volume v exc i in the isotropic phase consists in 
performing N att attempts of inserting two chains of length I in a box of volume V as described 
in the following: 

1. Set the counter N ov = 

2. Build first chain of length / randomly, according to the following procedure: 

(a) Insert a first randomly oriented monomer. 

(b) Insert a monomer M. bonded to a free site S on chain ends (S can be chosen 
randomly among the two free sites of the partial chain). The orientation of M 
will be random and its position will be chosen randomly within the available 
bonding volume between M. and S. The bonding volume between Ai and S is 
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defined as the volume corresponding to all possible center of mass positions of 
M with M bonded to S. 

(c) If the number of monomer inserted is I terminate otherwise go to 1). 

where the first monomer inserted is placed in the center of the box and it is oriented 
with its attractive sites parallel to the x-axis. 

3. Build a second chain of length /, where the first monomer inserted is placed randomly 
within the simulation box with a random orientation. 

4. Increase N ov by 1 if two monomers belonging to different chains overlap and the two 
chains are either not self-overlapping or forming a closed loop. 

5. if the number of attempts is less than N att go to 2) otherwise terminate. 
Then v exc i can be calculated as follows: 

v exd = (54) 

Natt 

A reasonable choice for the total number of attempts is N att = 10 6 . In a similar fashion 
one can also calculate the bonding volume [75] between two monomers. In this case one 
monomer is kept fixed in the center of the simulation box and the other one is inserted with 
random position and orientation for a total of N att attempts. The bonding volume will be: 

V b ^V (55) 

where the factor 4 accounts for the fact that two particles can form 4 different possible 
bonds and Nb on d is the number of times that the two monomers were bonded after a ran- 
dom insertion. Finally with the same procedure used to calculate the excluded volume in 
the isotropic phase we can evaluate the excluded volume in the nematic phase. The only 
difference is that now monomers have to be inserted with an orientation extracted from the 



Onsager angular distribution defined in Eq. (28), so that the excluded volume depends also 



on the parameter a. Again if N ov is the number of times that two monomers belonging to 
different clusters overlap and N att is the total number of attempts then we have: 

v excl (l,l,a) = (56) 
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XI. APPENDIX C 



In this Appendix we explain how to calculate the parameters An (a), k^(a) and B^(a) 
of the nematic free energy functional. As a preliminary step we check that v exc i(l,l',a) for 
a fixed value of a is a second order polynomial of I and V as assumed in Eq. (|5]). In Fig. [9] 
(a) we plot v exc i{l, I, ol) as a function of I for different values of a and Xq, v exc i(l, I, a) can be 
well represented by a parabolic function, in agreement with Eq. (j5j). 

We start by observing that the a dependence of A N (a), k N (a) and B N (a) in the case of 
hard cylinders following the Onsager distribution can be expanded in powers of a -1 / 2 as 



A N (a) 
k N (a) 
B N (a) 



c oo + 

ClO + 

C20 + 



Coi 
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C02 
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0)3 
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C04 


aVa 
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a 3 / 2 


o: 2 


Cll 
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C12 
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Cl3 
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C14 
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a 3 / 2 


a 2 


c 2 i 
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C22 
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C23 
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C24 


«V2 


a 


a 3 / 2 


a 2 



(57) 



where are the elements of the 3x4 matrix C. In the case of cylinders, some of the c^- 
vanishes|45j. We assume here that the same a dependence holds for SQ. 

In view of this result the co-volume as function of I and a can be expressed as 

du 
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UH) I v A J ^1 , ^2 . 



+ 



a 



(58) 



where c?^ p , for p = 0, 4 are fitting parameters. Fig|9] (b)-(d) shows the numerical calculation 
of the covolume varying a for three particular elongations (Xq = 1,2,3), together with fits 



to the functional form of Eq. (58). 

The good quality of the fits (reduced \ 2 is always much less than 1 for all fits) suggests 
that retaining terms up to 0(l/a 2 ) is to the present level of accuracy of our calculations 
absolutely appropriate. 

From these fits we can estimate the matrix C needed to evaluate the free energy in the 
nematic phase for each Xq. If we define in fact the following matrix P and the vectors q p , 
with p = ... 4 as follows: 



1 l a <- a 
1 h ll 
1 lr ll 



di bP 

\ d laP ) 



(59) 
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where l a , lb and l c are three different chain lengths for which we calculated the v exc i as a 
function of a, then we can calculate the matrix elements of C in the following way: 



( c \ 

VdClp 

\X$c 2p J 



p-v 



(60) 
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FIG. 9: (a) Excluded volume of two chains of length / as a function of chain length for the nematic 
cases a = 10,20,30,40 and three different elongations Xq = 1,2,3. (b)-(d) Excluded volume in 
the nematic phase calculated numerically as a function of a for two chains of equal length I, where 
I = 2, 3, 4, composed of monomers with Xq = 1, 2, 3. 
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